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Abstract 

This paper deals with the construction of a metamodel (i.e. a simplified mathematical 
model) for a stochastic computer code (also called stochastic numerical model or stochastic 
simulator), where stochastic means that the code maps the realization of a random variable. 
The goal is to get, for a given model input, the main information about the output probability 
distribution by using this metamodel and without running the computer code. In practical 
applications, such a metamodel enables one to have estimations of every possible random variable 
properties, such as the expectation, the probability of exceeding a threshold or any quantile. 
The present work is concentrated on the emulation of the quantile function of the stochastic 
simulator by interpolating well chosen basis function and metamodeling their coefficients (using 
the Gaussian process metamodel). This quantile function metamodel is then used to treat a 
simple optimization strategy maintenance problem using a stochastic code, in order to optimize 
the quantile of an economic indicator. Using the Gaussian process framework, an adaptive design 
method (called QFEI) is defined by extending in our case the well known EGO algorithm. This 
allows to obtain an “optimal” solution using a small number of simulator runs. 

Keywords: Adaptive Design, Asset management, Gomputer experiments. Expected Improve¬ 
ment, Gaussian process. Optimization, Stochastic simulator. Uncertainty. 


1 Introduction 

EDF looks for assessing and optimizing its strategic investments decisions for its electricity 
production assets by using probabilistic and optimization methods of “cost of maintenance 
strategies”. In order to quantify the technical and economic impact of a candidate maintenance 
strategy, some economic indicators are evaluated by Monte Garlo simulations using the VME 
software developed by EDF R&D (Lonchamp and Fessart jH]). The major output result of the 
Monte Garlo simulation process from VME is the probability distribution of the Net Present 
Value (NPV) associated to the maintenance strategy. From this distribution, some indica¬ 
tors, such as the NPV mean, the NPV standard deviation or the regret investment probability 
{¥{NPV < 0)) can easily be derived. 

Once these indicators have been obtained, one is interested in optimizing the strategy, for 
instance by determining the optimal investments dates leading to the highest mean NPV and the 
lowest regret investment probability. Due to the discrete nature of the events to be optimized, the 
optimization method is actually based on genetic algorithms. However, during the optimization 
process, one of the main issues is the computational cost of the stochastic simulator to optimize, 
which leads to methods requiring a minimal number of simulator runs (Dellino and Melon! [5]). 
Genetic algorithms require often several thousands of simulator runs and, in some cases, are no 
more applicable. 

The solution investigated in this study is a first attempt to break the computational cost of 
this problem by the way of using a metamodel instead of the simulator within the mathematical 
optimization algorithm. From a first set of simulator runs (called the learning sample and 
coming from a specific design of experiments), a metamodel consists in approximating the 
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simulator outputs by a mathematical model (Fang et al. [7]). This metamodel can then be 
used to predict the simulator outputs for other input configurations. 

Many metamodeling techniques are available in the computer experiments literature (Fang 
et al. [7]). Formally, the function G representing the computer model is defined as 

G: F; ^ 

X !->■ G{x) ' ' 

where E C (d € N*) is the space of the input variables of the computer model. Metamodeling 
consists in the construction of a statistical estimator G from an initial sample of G values 
corresponding to a learning set x with xG E and 4fx < oo (with the cardinal operator). 

However, the conventional methods are not suitable in the present framework because of the 
stochastic nature of the simulator: the output of interest is not a single scalar variable but a 
full probability density function (or a cumulative distribution function, or a quantile function). 
The computer code G is therefore stochastic: 

G : Exn ^ R , . 

(x,w) I—>■ G{x^u)) 

where H denotes the probability space. In the framework of robust design, robust optimization 
and sensitivity analysis, previous works with stochastic simulator consist mainly in approximat¬ 
ing the mean and variance of the stochastic output (Bursztyn and Steinberg [3, Kleijnen [S], 
Ankenman et al [T], Marrel et al. [ig, Dellino and Melon! i)- Other specific characteristics of 
the distribution of the random variable G(x) (as a quantile) can also be modeled to obtain, in 
some cases, more efficient algorithms (Picheny et al. [IH]). 

In this paper, as first proposed by Reich et al. [IS] (who used a simple Gaussian mixture 
model), we are interested in a metamodel which predicts the full distribution of the random 
variable G(x), Vx S E. We then define the following deterministic function /: 


f : E ^ E 

X ^ fx density of the r.v. G(x) 
with E the family of probability density functions: 


E = 


S F^(IR), positive, measurable. 



(3) 

(4) 


For a point x G E, building fx with the kernel method requires a large number realiza¬ 

tions of G(x). A recent work (Moutoussamy et al. (TS]) has proposed kernel-regression based 
algorithms to build an estimator / of the output densities, under the constraint that each call 
to / is cpu-time expensive. Their result stands as the starting point for the work presented in 
this paper (based on a Master internship report, see Browne jl]). 

The next section briefly presents the VME application case. In the third section, we propose 
to use the Gaussian process metamodel and develop an algorithm to emulate the quantile func¬ 
tion instead of the probability density function. In the fourth section, this metamodel is used 
to treat the VME case, in order to optimize a quantile. Using the Gaussian process metamodel 
framework, an adaptive design method is also proposed, allowing to obtain an “optimal” solu¬ 
tion using a small number of VME simulator runs. A conclusion synthesizes the main results of 
this paper. 


2 The VME application case 

On an industrial facilities, we are interested by the replacement cost of four components in 
function of the date of their maintenance (in year) (see Lonchamp and Fessart |llj for more 
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details). To these four inputs, one additional input is the date (in year) of recovering a new 
component. This input space is denoted F : 


F= (^{41,42,...,50} X {11,12,...,20}. 


F is therefore a discrete set (#F = 10® ). We have 


(5) 


G : 

F 


K 


X = (a:i,a;2,a;3,a;4,X5) 


NPV(a:), 

f ■■ 

F 

-> 

F 


X = {xi,X 2 ,X 3 ,X 4 ,X 5 ) (density of NPV(a:)). 

Figure [2 provides examples of the output density of VME. The 10 input values inside F have 
been randomly chosen. It shows that there is a small variability between the curves. 



Supp 


Figure 1: 10 ouput probability densities of the VME code (randomly sampled). 


The optimization process consists in finding the point of F which gives the “maximal NPV” 
value. As NPV(a;) is a random variable, we have to summarize its distribution by a deterministic 
operator FI, for example: 

H{g) = E(g) Wg G F (7) 

or 

H{g) = qg{a) Vg S JF (8) 

with qg(ot) the a-quantile of g. Our VME-optimization problem turns then to the determination 
of 

X* ;= arg max H{fV)- (9) 

xeF 

However, several difficulties occur: 

• VME is a cpu-time costly simulator and the size of the set F is large. Computing {fx)x^ f, 
needing NyiQ x ^F simulator calls (where A^mc is worth several thousands), is therefore 
impossible. Our solution is to restrict the VME calls to a learning set x C F (with 
#X = 200 in our application case), randomly chosen inside F. We will then have {fx)xex- 
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• Our metamodel, that we will denote 


{I 


xGF 


, cannot be computed on F due to its large 


size. A new set E is therefore defined, with E C F and = 2000. We will then have 
In this work, we limit our study to this restricted space E instead of the full 


(A) 


x£E 


space E. Other work will be performed to extend our methodology to the E study. 


3 Gaussian process metamodel of a stochastic simulator 

3.1 Basics on the Gaussian process model 

Let us consider n evaluations of a deterministic computer code. Each evaluation Y{x) G K 
of a simulator scalar output comes from a d-dimensional input vector x = (xi,... ,Xd) G E, 
where if is a bounded domain of The n points corresponding to the code runs are called 
the experimental design and are denoted as Xg = ... ,x^^'>). This learning sample comes 

from the learning set, previously noted x (n = #x)- Th® outputs will be denoted as Ys = 

■ ■ ■, with 2 /^*) = G'(x^®)) V i = 1. .n. Gaussian process modeling (Sacks et al. [20]), 
also called kriging model (Stein [21]), treats the simulator deterministic response G(x) as a 
realization of a random function Y(x), including a regression part and a centered stochastic 
process: 

Y (x) = h{x) + Z (x). (10) 

The deterministic function h{x) provides the mean approximation of the computer code. We 
can use for example a one-degree polynomial model: 

d 

h{x) = ^o + '^l3jXj , (11) 

i=i 

where /3 = [/3o, • ■ ■, PdY is the regression parameter vector. The stochastic part Z{x) is a Gaus¬ 
sian centered stationary process fully characterized by its covariance function: Gov(Z(x), Z(u)) = 
a^Ke{x — u), where denotes the variance of Z, Kg is the correlation function and 9 gM.^ is 
the vector of correlation hyperparameters. This structure allows to provide interpolation and 
spatial correlation properties. Several parametric families of correlation functions can be chosen 
(Stein [2T]L 

If a new point x* = (x*,... ,x2) G E is considered, we obtain the predictor and variance 
formulas for the scalar output Y{x*): 

E[r(x*)|n] = h{x*) + fc(x*)*s;i(n - h{x,)) , 

MSE{x*) = Var[y(x*)|y«] = - k{x*YT,J^k{x*) , 

with 

k{x*) = [Gov( 2 /(i),y(x*)),...,Gov( 2 /(-),y(x*))]‘ 

= aYKe{x(^\x*),...,Ke{x^^\x*))Y 

and the covariance matrix 

The conditional mean (Eq. @) is used as a predictor. The variance formula (Eq. (HI) 
corresponds to the mean squared error (MSE) of this predictor and is also known as the kriging 
variance. This analytical formula for MSE gives a local indicator of the prediction accuracy. 
More generally, Gaussian process model defines a Gaussian distribution for the output variable 
at any arbitrary new point. This distribution formula can be used for uncertainty and sensitivity 
analysis (Marrel et al. [TT], Le Gratiet et al. (TD])- Regression and correlation parameters /3, a 
and 9 are ordinarily estimated by maximizing likelihood functions (Fang et al. [7]). 


( 12 ) 

(13) 

(14) 

(15) 
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3.2 Emulation of the simulator quantile function 
3.2.1 General principles 

In our VME-optimization problem, we are especially interested by several quantiles (for exam¬ 
ple at the order 1%, 5%, 50%, 95%, 99%) rather than statistical moments. In Moutoussamy 
et al. m and Browne [1], quantile prediction with density-based emulator has shown some 
deficiencies. Therefore, instead of studying Eq. we turn our modeling problem to 

G : E ^ M 

X = {xi,X2,X3,X4,X5) NPV(a;), 

(16) 

Q : E —>■ Q 

X ^ Qx quantile function of NPV(a;) 

where Q is the space of increasing functions defined on ]0,1[, with values in [a, b] (which is the 
support of the NPV output). For x G E, a quantile function is defined by : 

VpG]0,l[, Qxip) = t G [a,b] such as f fxis)de = p. (17) 

J a 

For the same points than in Figure Figure shows the 10 quantile function ouputs Q which 
present a rather low variability. 



SuppQ[ 6 : 96 ] 


Figure 2: Quantile functions Q for 10 points of E (randomly chosen). 


We consider the learning set x (n = =ffx) A^mc ^ xl G-simulator calls in order to obtain 
, the empirical quantile functions of (NPV(a;)),j,g^. In this work, we will use IVmc = 


10^, which is sufficiently large to obtain a precise estimator of Q^c with Qx^'^ ■ Therefore, we 
neglect this Monte Carlo error. In the following, we simplify the notations by replacing Qx^'^ 
by Qx- 

The approach we adopt is similar to the one used in metamodeling a functional output of a 
deterministic simulator (Bayarri et ah, [2], Marrel et. al. m)- The first step consists in finding 
a low-dimensional functional basis in order to reduce the output dimension by projection, while 
the second step consists in emulating the coefficients of the basis functions. However, in our 
case, due to the nature of the functional outputs (quantile functions), some particularities will 
arise. 
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3.2.2 Projection of Qx by the Modified Magic Points (MMP) algorithm 

Adapted from the Magic Points algorithm (Maday et al. [T7j) for probability density functions, 
the MMP algorithm has been proposed in Moutoussamy et al. m- It is a greedy algorithm 
that builds interpolator (as a linear combination of basis functions) for a set of functions by 
iteratively picking a basis function in the learning sample output set and a set of interpolation 
points. At each step j G {2 ,... ,9} of the construction of the functional basis, one picks the 
element of the learning sample output set that maximizes the gap (L^ distance) between this 
element and the interpolator which used the previous j — 1 functions of the basis. The total 
number q of functions is chosen with respect to a convergence criterion. Mathematical details 
will not be provided in the present paper. 

In this paper, we apply the MMP algorithm on quantile functions. The first step consists in 
a projection of {Qx)xex' 

Qx = '^’tpjix)Rj'ix G X (18) 

i=i 

where '0 = {fpi, ■ ■ ■ ,fpq) (the coefficients) and R = {Ri,...,Rq) (the quantile functions of the 
basis) are determined by MMP. We need to restrict the solutions to the following constrained 
space: 

C={0eK'J, 01,...,0, >0} (19) 

in order to ensure the monotonic increase of Qx- 

In our VME application, the choice q = 5 has shown sufficient approximation capabilities. 
For one example of quantile function output, a small relative L^-error (0.2%) between the 
observed quantile function and the projected quantile function is obtained. Figure confirms 
also the relevance of the MMP method. 



SuppQ 

Figure 3: For one point x G Xi Qx (red line) and Qx (black points). 


3.2.3 Gaussian process metamodeling of the basis coefficients 

Estimations of the 0(a;) = (0i(x),... ,0g(x)) {x G E) coefficient vector will be performed with 
q Gaussian process metamodels in order to build Qx'- 

q 

Qx = Y^ '^jix)Rj 'ix G E- (20) 


6 




To ensure that ipj £ C (j = 1... g), we use a logarithmic transformation: 


Ti : C 

tjj 


and its inverse transformation: 

Ta : 


{l0g{tpl + 1), ...,l0g(V'g + 1)) 


c 

(exp((;ii) - 1, - 1). 


( 21 ) 


( 22 ) 


We then compute (/)(x) := Ti{ip{x)) \/x G x and suppose that 4> is a Gaussian process realization 
with q independent margins, (j) is estimated by 


^{x) := 'E\Yx I Ys] \/x G E 


(23) 


with Yg the learning sample output. We obtain 

'tpix) := T2{^{x)) \/x G E (24) 


and Eq. (20) can be applied as our metamodel predictor of the quantile function. 

In our VME application, we have built the metamodel on the set E (with the choice q = 5). 
For one example of quantile function output, a small relative L^-error (2.8%) between the 
observed quantile function and the emulated quantile function is obtained. Figure confirms 
also the relevance of the metamodeling method. 



Figure 4: For one point x G Xi Qx (red line) and Qx (black points). 


4 Application to an optimization problem 


4.1 Direct optimization on the metamodel 

We now apply our quantile function metamodel with a quantile-based objective function 


H : Q ^ R 
q q{p) 


with p s]0,1[. We look for 


X* := arg max Qa; (p) 
xeF 


(25) 

(26) 
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but have only access to 


(27) 


X* := arg max Qa; (p). 

x£F 

We study the relative error of H{Q) on E by computing 

. XX ^ ^ X XX X ( X! I ‘5x(p) - Qx{p) I I ■ (28) 

Taiax^(zE{Qx{p))-^-^^xGEiQxip)) ) 

As an example, for p = 0.5 (median estimation), we find 

max^(zE{Qx{p)) = 0.82, maxx(zE{Qx{p)) = 0.42, err = 5.4%. (29) 

If we define y = argmax^.^^ Qx{p) the best point from the metamodel, we obtain Qy{p) = 0.29 
while vaaxx^xQxip) = 0.35. The exploration of E by our metamodel does not bring any 
information. We have observed the same result by repeating 100 times the experiments (changing 
the initial design). It means that the punctual errors on the quantile function metamodel are 
too large for this optimization algorithm. In fact, the basis functions i?i, ...,i ?5 that the MMP 
algorithm has chosen on y are not able to represent the extreme parts of the quantile functions 
of A. 

As a conclusion of these tests, the quantile function metamodel cannot be directly applied 
to solve the optimization problem. In the next section, we propose an adaptive algorithm 
which consists in sequentially adding simulation points in order to capture interesting quantile 
functions to be added in our functional basis. 


4.2 QFEI: An adaptive optimization algorithm 

After the choice of x, E and the families (Qa;)^,^^, (Ox)xgx (Ox)xgBi our new algorithm 
will propose to perform new interesting (for our specific problem) calls to the VME simulator 
on E (outside of x). With the Gaussian process metamodel, which provides a predictor and 
its uncertainty bounds, this is a classical approach used for example in black-box optimization 
problem (Jones et al. [5]) and rare event estimation (Beet et al. [3]). The goal is to provide 
some algorithms which mix global space exploration and local optimization. 

Our algorithm is based on the so-called EGO (Efficient Global Optimization) algorithm 
(Jones et al. [S]) which uses the Expected Improvement (El) criterion to optimize a deterministic 
simulator. Our case is different as we want to maximize 


H : E ^ R 

X !--)■ Qx{p) p-quantile of NPV(x). 


(30) 


We will then propose a new algorithm called the QFEI (for Quantile Function Expected Im¬ 
provement) algorithm. 

As previously, we use the set E C F with ^E = 5000 {E is a random sample in F), the 
initial learning set x C F with ^x = 200 (initial design of experiment), {Qx)x^x’ 

and {Qx)x^e- We denote D the current learning set (the initial learning set increased with 
additional points coming from QFEI). As Gaussianity will be needed on the components of ip, 
we did not performed a logarithmic transformation as in Section |3.2.3| In our case, it has not 
implied negative consequences. 

We apply the Gaussian process metamodeling on the q independent components ipi, ...jipq- 


ipj (x) ~ M 


>j{x),MSEj{x) j Vj G {1,..., q} Va; G E. 


(31) 





Va; G E. 


(32) 


Qx{p) ^ M \^x{p),y^^Rj{pYMSEj{x) 

Then Qx{p) is a realization of the underlying Gaussian process Ux = J2‘j=i'^ji^)Rjip) with 


Ud 

{Ux)xeD 
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Ux 

:= E[Ux 1 Ud] 

Vx G E, 

(33) 


:= Yar[Ux \ Ud] 

Vx G E. 



The conditional mean and variance of Ux are directly obtained from the q Gaussian process 
metamodels of the ip coefficients. 

At present, we propose to use the following improvement random function: 

I : E ^ R 

X i-t (C/a: — max (Gd))"'’. 

In our adaptive design, finding a new point consists in solving: 

Xnew := argmaxE[/(x)]. 
xGE 

Added points are those which have more chance to improve the current optimum. The expec¬ 
tation of the improvement function writes (the simple proof is given in Browne |3]): 

E[J(a:)] = (Ti/\d(x) {u{x)(j){u{x)) + ip{u{x))) \/x G E, with u{x) = — (3®) 

<xu\d{x) 

where p and (j) correspond respectively to the density and distribution functions of the reduced 
centered Gaussian law. 

In practice, several iterations of this algorithm are performed, allowing to complete the 
experimental design D. At each iteration, a new projection functional basis is computed and the 
q Gaussian process metamodels are re-estimated. The stopping criterion of the QFEI algorithm 
can be a maximal number of iterations or a stabilization criterion on the obtained solutions. No 
garantee on convergence of the algorithm can be given. In conclusion, this algorithm provides 
the following estimation of the optimal point x*: 


(34) 

(35) 


X* := argmax([/D), 

x^D 


(37) 


In our application case, we have performed all the simulations in order to know {Qx)x£Ei 
therefore the solution x*. Our first objective is to test our proposed algorithm for p = 0.4 which 
has the following solution: 


X* =(41,47,48,45,18) 
Q,.(p) =-1.72. 


We have also computed 
1 


— ^ Qx{p) = -3.15, Var {{Qx^x^e) = 0-59. 

^ xGE 


(38) 


(39) 


We start with D := x and we obtain 


max (Qx) = —1.95 


(40) 




(41) 


After 50 iterations of the QFEI algorithm, we obtain: 

f X* =(41,47,45,46,19) 

1 Q**(p) =-1.74. 

We observe that x* ~ x* and (p) — Qx* {p) which is a first confirmation of the relevance 
of our method. With respect to the initial design solution, the QFEI has allowed to obtain a 
strong improvement of the proposed solution. 50 repetitions of this experiment (changing the 
initial design) has also proved the robustness of QFEI. The obtained solution is always one of 
the five best points on E. 

QFEI algorithm seems promising but a lot of tests remain to perform and will be pursued 
in future works: changing p (in particular testing extremal cases), increasing the size of E, 
increasing the dimension d of the inputs, ... 

5 Conclusion 

In this paper, we have proposed to build a metamodel of a stochastic simulator using the 
following key points: 

1. Emulation of the quantile function which proves better efficiency for our problem than the 
emulation of the probability density function; 

2. Decomposition of the quantile function in a sum of the quantile functions coming from the 
learning sample outputs; 

3. Selection of the most representative quantile functions of this decomposition using an 
adaptive choice algorithm (called the MMP algorithm) in order to have a small number of 
terms in the decomposition; 

4. Emulation of each coefficient of this decomposition by a Gaussian process metamodel, by 
taking into account constraints ensuring that a quantile function is built. 

The metamodel is then used to treat a simple optimization strategy maintenance problem 
using a stochastic simulator (VME), in order to optimize an output (NPV) quantile. Using the 
Gaussian process metamodel framework and extending the El criterion to quantile function, the 
adaptive QFEI algorithm has been proposed. In our example, it allows to obtain an “optimal” 
solution using a small number of VME simulator runs. 

This work is just a first attempt and needs to be continued in several directions: 

• Gonsideration of a variable whose decrease could help to fight against the computa¬ 
tional cost of the stochastic simulator, 

• Improvement of the initial learning sample choice by replacing the random sample by a 
space filling design (Fang et al. [7j), 

• Algorithmic improvements to counter the cost of the metamodel evaluations and to increase 
the size of the study set E, 

• Multi-objective optimization (several quantiles to be optimized) in order to take advantage 
of our powerful quantile function emulator, 

• Application to more complex real cases, 

• Gonsideration of a robust optimization problem where environmental input variables of 
the simulator has not to be optimized but just create an additional uncertainty on the 
output. 
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